##########################################################################################
## https://sunduanchen.github.io/Scissor/vignettes/Scissor_Tutorial.html
library(data.table)
library(optparse)
library(ggplot2)
library(dplyr)
library(patchwork)
library(Seurat)
library(ggrepel)
library(parallel)

##########################################################################################

option_list <- list(
    make_option(c("--single_cell_file"), type = "character"),
    make_option(c("--cds_file"), type = "character"),
    make_option(c("--pathway_path"), type = "character"),
    make_option(c("--gene"), type = "character"),
    make_option(c("--singleCell_sample_file"), type = "character"),
    make_option(c("--single_cell_all_file"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){

    gene <- "MUC6"

    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic"
    single_cell_file <- paste0(work_dir,"/images/singleCell_MUC6/Scissor_STAD_MUC6_mutation.IM.CellRate.all.RData")
    out_path <- paste0("~/20220915_gastric_multiple/dna_combinePublic/images/singleCell_MUC6/diff_new")
    cds_file <- "~/ref/PCAWG_Elements/web_hg19/gc19_pc.cds.use.bed"
    pathway_path <- "~/ref/Pathway/"

    singleCell_sample_file <- "~/20220915_gastric_multiple/dna_combinePublic/config/singleCell_Sample.useThree.list"
    single_cell_all_file <- paste0(work_dir,"/public_ref/singleCell/njmu/epiall_nor_PCA_50_RE0.5.Rdata")

}

##########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

out_path <- opt$out_path
gene <- opt$gene
single_cell_file <- opt$single_cell_file
singleCell_sample_file <- opt$singleCell_sample_file
single_cell_all_file <- opt$single_cell_all_file
pathway_path <- opt$pathway_path
cds_file <- opt$cds_file

dir.create(out_path , recursive = T)

##########################################################################################
## 提取用到的4个肠化样本
info_singlecell <- data.frame(fread(singleCell_sample_file))

sc_dataset <- load(single_cell_all_file, verbose = F)
sc_dataset_all <- epiall_nor_PCA_50_RE0.5
##Idents函数定义要取的类别
Idents(sc_dataset_all) <- "sample"   

## 提取用到的样本
sc_patient_id <- unique(substring(info_singlecell$singlecell_ID ,1 ,5 ))
sc_dataset_im <- subset(sc_dataset_all , patient %in% sc_patient_id & sample=="IM" )
Idents(sc_dataset_im) <- "celltype" 
## 取pit
sc_dataset_im <- subset(sc_dataset_im , idents=c("Pit"))


##########################################################################################
## MUC6突变的样本
sc_datasets <- load(single_cell_file, verbose = F)
coding_gene <- data.frame(fread(cds_file))
coding_gene <- unique(sapply(strsplit(coding_gene$V4 , "::") , "[" , 3))
## sc_dataset
##Idents函数定义要取的类别
##########################################################################################
Idents(sc_dataset) <- "celltype" 
## 取pit
sc_dataset <- subset(sc_dataset , idents=c("Pit"))

## 比较scissor + 和其它细胞的差异
##   Scissor_select[infos4$Scissor_pos] <- 1
##   Scissor_select[infos4$Scissor_neg] <- 2
sc_dataset$scissor <- ifelse( sc_dataset$scissor == 1 , "Positive" , "Other" )
sc_dataset$scissor <- factor( sc_dataset$scissor , levels = c("Positive" , "Other") , order = T )

##########################################################################################
## 注释整体里面的突变情况
sc_dataset_im$MUC6_mut <- ifelse( sc_dataset_im$patient == "JZ732" , "Mut" , "Wild" )
positive_cell <- names(sc_dataset$scissor)[sc_dataset$scissor=="Positive"]
sc_dataset_im$scissor <- "Other"
sc_dataset_im$scissor[which(names(sc_dataset_im$celltype) %in% positive_cell)] <- "Positive"

##########################################################################################
## 差异表达,scissor阳性的Pit
Idents(sc_dataset_im) <- "scissor" 

## 不限制foldchange
## 显示所有基因
pit_deg <- FindMarkers(sc_dataset_im,
        ident.1 = "Positive" , ident.2 = "Other" , only.pos=F, min.pct=0, logfc.threshold=0
        )

## 提取coding基因
pit_deg_coding <- pit_deg[rownames(pit_deg) %in% coding_gene,]
pit_deg_coding$p_val_adj <- p.adjust( pit_deg_coding$p_val , method = "fdr" )

## 提取显著差异的基因
use_dat <- pit_deg_coding

#对原数据进行处理
use_dat$log10_q_Value <- -log10(use_dat$p_val_adj)
use_dat$log_FC <- use_dat$avg_log2FC
use_dat$gene <- rownames(use_dat)

#设置阈值
q_t <- 0.05
foldchange_t <- 1.5
logFC_cutoff <- log2(foldchange_t)
log10_P_Value_cutoff <- -log10(q_t)

data <- use_dat[,c('gene','log_FC','log10_q_Value')]
colnames(data) <- c('gene','log_FC','-log10_P_Value')
image_name <- paste0( out_path , "/DiffGene.Pit.scissorP_vs_OtherPit.tsv" )     
write.table( data , image_name , row.names = F , sep = "\t" , quote = F )

## Mut中显著高表达的基因
image_name <- paste0( out_path , "/DiffGene.Pit.high.scissorP_vs_OtherPit.tsv" )     
write.table( subset( use_dat , log_FC >= logFC_cutoff & log10_q_Value > log10_P_Value_cutoff ) , image_name , row.names = F , sep = "\t" , quote = F )

## 展示想展示的基因
data$gene <- ifelse( data$gene %in% c("GKN1" , "GKN2" , "DSC2" , "DSG2") , data$gene , "" )

col <- c('#376B6D','#9F353A' , '#F7C242','#BDC0BA' )
names(col) <- c('DOWN', 'UP', 'TARGET', "GREY" )

alpha <- 1

options(ggrepel.max.overlaps = 200)
plot1 <- ggplot(data = data,aes(x = log_FC,y = `-log10_P_Value`))+
        geom_point(data = subset(data,abs(log_FC)<logFC_cutoff),
                aes(size = abs(log_FC)),col = col["GREY"],alpha = alpha)+
        geom_point(data = subset(data,abs(`-log10_P_Value`)<log10_P_Value_cutoff & abs(log_FC)>logFC_cutoff),
                aes(size = abs(log_FC)),col = col["GREY"],alpha = alpha)+
        geom_point(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & log_FC>logFC_cutoff),
                aes(size = abs(log_FC)),col = col["UP"],alpha = alpha)+
        geom_point(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & log_FC< -logFC_cutoff),
                aes(size = abs(log_FC)),col = col["DOWN"],alpha = alpha)+
        geom_point(data = subset(data,gene!=""),
                aes(size = abs(log_FC)),col = col["TARGET"],alpha = alpha)+
        theme_bw()+
        labs(x='log fold change\n(Scissor+ cells/All other cells in three samples )',y='-log10(adjusted p-value)')+
        geom_vline(xintercept = c(-logFC_cutoff,logFC_cutoff),lty = 3,col = 'black',lwd = 0.4)+
        geom_hline(yintercept = log10_P_Value_cutoff,lty = 3,col = 'black',lwd = 0.4) +
        geom_text_repel(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & abs(log_FC)>logFC_cutoff), 
            fontface="bold", nudge_x = 0.65, nudge_y = 1, 
            aes(label = gene),size = 5,col = 'black' , face = 'bold') +
        theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='right',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                legend.text = element_text(size = 12,color="black",face='bold'),
                axis.text.x = element_text(size = 12,color="black",face='bold'),
                axis.text.y = element_text(size = 10,color="black",face='bold'),
                axis.title.x = element_text(size = 14,color="black",face='bold'),
                axis.title.y = element_text(size = 14,color="black",face='bold'),
                strip.text.x = element_text(size = 12,color="black",face='bold'),
                axis.ticks.length = unit(0.2, "cm") ,
                axis.line = element_line(size = 0.5))

image_name <- paste0( out_path , "/DiffGene.Pit.scissorP_vs_OtherPit.pdf" )     
ggsave( image_name , plot1 , width = 6 )


##########################################################################################
## 差异表达，MUC6突变患者的Pit
Idents(sc_dataset_im) <- "MUC6_mut" 

## 不限制foldchange
## 显示所有基因
pit_deg <- FindMarkers(sc_dataset_im,
        ident.1 = "Mut" , ident.2 = "Wild" , only.pos=F, min.pct=0, logfc.threshold=0
        )

## 提取coding基因
pit_deg_coding <- pit_deg[rownames(pit_deg) %in% coding_gene,]
pit_deg_coding$p_val_adj <- p.adjust( pit_deg_coding$p_val , method = "fdr" )

## 提取显著差异的基因
use_dat <- pit_deg_coding

#对原数据进行处理
use_dat$log10_q_Value <- -log10(use_dat$p_val_adj)
use_dat$log_FC <- use_dat$avg_log2FC
use_dat$gene <- rownames(use_dat)

#设置阈值
q_t <- 0.05
foldchange_t <- 1.5
logFC_cutoff <- log2(foldchange_t)
log10_P_Value_cutoff <- -log10(q_t)

data <- use_dat[,c('gene','log_FC','log10_q_Value')]
colnames(data) <- c('gene','log_FC','-log10_P_Value')
image_name <- paste0( out_path , "/DiffGene.Pit.MUC6Mut_vs_OtherThree.tsv" )     
write.table( data , image_name , row.names = F , sep = "\t" , quote = F )

## Mut中显著高表达的基因
image_name <- paste0( out_path , "/DiffGene.Pit.high.MUC6Mut_vs_OtherThree.tsv" )     
write.table( subset( use_dat , log_FC >= logFC_cutoff & log10_q_Value > log10_P_Value_cutoff ) , image_name , row.names = F , sep = "\t" , quote = F )

## 展示想展示的基因
data$gene <- ifelse( data$gene %in% c("GKN1" , "GKN2" , "DSC2" , "DSG2") , data$gene , "" )

col <- c('#376B6D','#9F353A' , '#F7C242','#BDC0BA' )
names(col) <- c('DOWN', 'UP', 'TARGET', "GREY" )

alpha <- 1

options(ggrepel.max.overlaps = 200)
plot1 <- ggplot(data = data,aes(x = log_FC,y = `-log10_P_Value`))+
        geom_point(data = subset(data,abs(log_FC)<logFC_cutoff),
                aes(size = abs(log_FC)),col = col["GREY"],alpha = alpha)+
        geom_point(data = subset(data,abs(`-log10_P_Value`)<log10_P_Value_cutoff & abs(log_FC)>logFC_cutoff),
                aes(size = abs(log_FC)),col = col["GREY"],alpha = alpha)+
        geom_point(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & log_FC>logFC_cutoff),
                aes(size = abs(log_FC)),col = col["UP"],alpha = alpha)+
        geom_point(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & log_FC< -logFC_cutoff),
                aes(size = abs(log_FC)),col = col["DOWN"],alpha = alpha)+
        geom_point(data = subset(data,gene!=""),
                aes(size = abs(log_FC)),col = col["TARGET"],alpha = alpha)+
        theme_bw()+
        labs(x='log fold change\n(Scissor+ cells/All other cells)',y='-log10(adjusted p-value)')+
        geom_vline(xintercept = c(-logFC_cutoff,logFC_cutoff),lty = 3,col = 'black',lwd = 0.4)+
        geom_hline(yintercept = log10_P_Value_cutoff,lty = 3,col = 'black',lwd = 0.4) +
        geom_text_repel(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & abs(log_FC)>logFC_cutoff), 
            fontface="bold", nudge_x = 0.65, nudge_y = 1, 
            aes(label = gene),size = 5,col = 'black' , face = 'bold') +
        theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='right',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                legend.text = element_text(size = 12,color="black",face='bold'),
                axis.text.x = element_text(size = 12,color="black",face='bold'),
                axis.text.y = element_text(size = 10,color="black",face='bold'),
                axis.title.x = element_text(size = 14,color="black",face='bold'),
                axis.title.y = element_text(size = 14,color="black",face='bold'),
                strip.text.x = element_text(size = 12,color="black",face='bold'),
                axis.ticks.length = unit(0.2, "cm") ,
                axis.line = element_line(size = 0.5))

image_name <- paste0( out_path , "/DiffGene.Pit.MUC6Mut_vs_OtherThree.pdf" )     
ggsave( image_name , plot1 , width = 6 )
